% 参数设置
global A B
G = 6.67430e-11; % 万有引力常量
mo = 1.989e30; % 太阳质量
c = 299792458; % 光速
h = 2 * 1.32712440018e20; % 水星掠面速度的两倍
% 初始条件
A=G*M/h^2;B=3*G*M/c;

[thi,u]=ode45(@sxjd,0:0.001*pi:30*pi,[0.2,0]);
[x,y]=pol2cart(thi,1./u(:,1));
comet(x,y,0.05); 
% 绘制轨道
plot(0,0,'k*');
plot(x,y,'k')
hold on
axis equal

function F=sxjd(thi,u)
global A B
F=[u(2);
    -u(1)+A+B*u(1)^2];
end
